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Abstract The purpose of this study is to evaluate the 
effect of growth rate on intra-tree variation in basic density 
of hinoki cypress (Chamaecyparis obtusa ) quantitatively 
using the statistical modeling technique. Nineteen sample 
trees were harvested from 50-year-old hinoki stand which 
consists of two different growth rate plots. Disks were cut 
from sample trees at height positions of 2, 4 m, and then 
4 m intervals until 16 m position. Radial strips were cut 
from the disks, and ring widths and basic density were 
measured at 5-ring intervals. The basic density decreased 
with age at any height positions. The linear mixed model 
was fitted to the age trend data having two nested grouping 
levels, i.e., tree and position within tree. Models having 
various mean and covariance structures were tested in 
devising an appropriate wood density model. The model, 
consisting of the mean structure with quadratic function of 
cambial age was able to describe the intra-tree variation in 
basic density. The model containing the random effects 
which consist of effect of the tree level and vertical stem 
position level explained the density variation adequately. 
The growth rate did not show the significant effect on the 
basic density variation within the stem. 
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Introduction 

Hinoki cypress (Chamaecyparis obtusa) is one of the most 
important commercial species in Japan and its afforestation 
area accounts for 25 % of all plantation forest area in Japan 
[1]. Historically, hinoki has been used as building materials 
because it shows remarkable mechanical properties and 
durability [2]. For instance, it is well known that Horyuji, 
the world’s oldest wooden architecture, was built more 
than 1,300 years ago and hinoki was mainly used as its 
construction member [3]. However, the wood resources 
have shifted from natural forests to plantation forests 
recently, and thus the wood properties of hinoki from 
plantation forests would be different from the past situa¬ 
tion. The correct recognition of the wood properties in 
current hinoki resources should be necessary for the ade¬ 
quate forest management and utilization. 

Numerous studies have reported the variation of wood 
properties in a number of species [4]. However, many of 
these reports are not quantitative examinations but quali¬ 
tative examinations, except for the genetic variation which 
provides various genetic parameters. If the variation pat¬ 
terns of wood properties can be formulated quantitatively, 
one can interpolate the data or extrapolate the future situ¬ 
ations based on the formulae. 

There may be no discussed matter with regard to wood 
quality variation than relationship between growth rate and 
wood properties [4]. Especially, the effects of growth rate 
on wood density have been studied intensively. In hinoki 
cypress, Hirai reported that high growth rate would 
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produce wood having low basic density [5]. Fujiwara et al. 
[6] examined the variation of basic density of hinoki which 
was obtained from nine test stands, and reported that there 
were no significant differences between thinned stands and 
un-thinned stands. Growth rate affects wood properties at 
different tree ages, and ignoring this effect could result in 
growth rate being incorrectly identified as the cause of 
differing wood properties [7, 8]. 

Sequences of yearly measurements of wood character¬ 
istics are considered as longitudinal (time series) data, and 
the analysis of longitudinal data is able to evaluate the age 
effect on variation of these characteristics. In this case, 
these data are repeated measures made of the same char¬ 
acteristic on the same observational unit (tree, disk, and 
ring), and such data generally present temporal autocorre¬ 
lation, heteroscedasticity, and nonstationarity of the mean 
[9-11]. The mixed-effects analysis technique is frequently 
used for such grouped data including longitudinal data, 
repeated measures data, and multilevel data [12, 13]. 
Correlations among observations made on the same subject 
or experimental unit can be modeled using random effects, 
random regression coefficients, and through the specifica¬ 
tion of a covariance structure. Recent studies have devel¬ 
oped models predicting variation in wood properties within 
a stem or cross section using the mixed-effects analysis for 
many species, such as pine and spruce [14-21]. 

Although the general variation pattern of basic density 
within stem is well known in hinoki, there is little infor¬ 
mation available about the effects of growth rate on the 
variation pattern of basic density. The objective of this 
study is to evaluate the effect of growth rate on variation 
pattern of basic density within stem. The linear mixed- 
effects model was used to explain the variation pattern of 
basic density quantitatively, and the effects of growth rate 
on its variation were also assessed on the basis of statistical 
manner. 

Materials and methods 

Sample materials were obtained from 50-year-old hinoki 
cypress stands in Tottori University Forest located in 
Maniwa, Okayama (35° 16' N, 133°36 / E; approximately 
540 m elevation). The annual average precipitation and 
temperature from 2002 to 2012 in the research forest were 
1989 mm and 11.6 °C, respectively. Hinoki cypress seed¬ 
lings were planted in 1962 at an initial plant density of 
3333 trees/ha. The test stand consisted of four 10-m square 
plots, which were two fast growth plots and two slow 
growth plots. They were adjacent to each other. The 
information of the test stand is summarized in Table 1. 
Although the detailed archive of the stand is unclear, the 
fast growth plots had been thinned several times. On the 


Table 1 Characteristics of the plots and sample trees in the test stand 



Plot A 


Plot B 


Al 

A2 

B1 

B2 

Stand density (tree/ha) 

1000 

1200 

2100 

2800 

BA (m 2 /ha) 

50.9 

52.5 

63.2 

42.0 

DBH (cm) 

22.8 

20.9 

19.9 

13.6 

Height (m) 

17.6 

17.5 

17.2 

11.5 

BLC (m) 

11.2 

9.1 

11.1 

5.5 


The data were obtained when the sample trees were harvested. DBH, 
Height, and BLC are the mean value of sample trees 

BA basal area, DBH diameter at breast height, BLC, height to the base 
of live crown 

other hand, the slow growth plots have never been thinned 
previously. As a result, the stand density between them was 
quite different. 

From each plot, four or five, totally nineteen sample 
trees, were felled in 2012. Height and height to the base of 
live crown (BLC) of the sample trees were measured. The 
average height, DBH and BLC of sample trees were 

15.9 m, 20.7 cm and 9.1 m, respectively. The height ran¬ 
ged from 10.2 to 18.6 m. DBH ranged from 11.9 to 

30.9 cm. Analysis of variance confirmed the significant 
difference of the DBH among plots (p = 0.006; data not 
shown). 

A 5-cm-thick, knot-free sample disks were cut from 
each sample tree at height positions of 2, 4 m, and then 4 m 
intervals until 16 m position. A 3-cm-thick radial strip was 
cut from each disk and ring width was measured from pith 
to outward in every five rings. Radial diameter of sample 
strips ranged from 1.4 to 12.1 cm. The measurements were 
carried out for two radial directions and ring width was 
expressed as the mean of both values. After that, the strips 
were cut into every 5 rings for basic density analysis. 
Water displace method was used to determine the basic 
density. The basic density was also determined as the mean 
of both radial directions. Radial and longitudinal variations 
of basic density for each plot are presented in Fig. 1. 

Model development 

The start model 

The mixed model technique [12, 13] was used for modeling 
the effects of growth rate on the intra-tree variation in basic 
density. Figure 1 indicates that basic density can be rep¬ 
resented by the quadratic function of cambial age, i.e., ring 
number from pith. On the start model, the fixed effects 
consisted of population mean and effects of plots deter¬ 
mined by growth rate. The random effects consisted of 
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Fig. 1 Radial and longitudinal variations of basic density from different growth rate plots 


effect of the tree level and vertical stem position level. The 
model expressed as 


BDhijk — (ft 0 + ft oh + boi + b 0lJ ) + (ft 1 + fiih + bu + bujj AGE\ 
+ (ft2 + ftlh + ^ 2 i + b2ij) AGE 2 + Ehijk , 



i 

^5- 

o 

_1 



1 

o 

i_ 

bi = 

1 

S3- ^5- 

to 

i_ 

* l), 

btj = 

i 

<n 

_i 


W 2 ), 


Shijk~N(0, a 2 ), 



where BD hijk was the basic density of the kth cambial age 
of the yth vertical stem position in the ith tree in the hth 
plot; ft 0 , fti and ft 2 were the population mean of the basic 
density; ft 0h , ft lh and ft 2 h were the fixed-effect parameters 
of hth plot; bi was the tree-level random-effects vector; b t j 
was the position-level random-effects vector; £ hi j k was the 
within-group error. The b t were assumed to be independent 
for different i, the b t j were assumed to be independent for 
different i, j and independent of the b h and the £ hijk were 
assumed to be independent for different i, j , k and inde¬ 
pendent of the random effects. The large number of 
parameters in Eq. (1) makes the optimization of the pro¬ 
filed log-restricted-likelihood quite difficult and unstable 
[12]. To make the optimization more stable during this 
model building phase, we simplify Eq. (1) by assuming *Pi 
and 2 as diagonal matrices. The models in this article 
were fitted using the nlme package in R version 3.0.0 [22]. 


Selecting the fixed-effects structure 

First, we evaluate whether the quadratic function of cam¬ 
bial age is adequate to describe the observed data and also 
test whether the growth rate has significant effect on the 
intra-tree variation on the basic density. The result of fitting 
indicated that both first- and second-order terms of age 
were highly significant (p < 0.001), and thus, the quadratic 
function of cambial age well describes the radial variation 
of basic density. 

There were no clear effects in the terms of plot 
(p = 0.087), plot-AGE interaction (p = 0.288), and plot- 
AGE interaction (p = 0.079). The results indicate that 
growth rate does not affect the variation pattern of basic 
density within the stem in hinoki cypress. Consequently, 
the fixed effect of the Eq. (1) could be simplified to the 
structure without plot effects. 

Determining the variance-covariance structure 
of random effects 

The age quadratic model without plot effects (model 1.1) 
was examined to determine the variance-covariance 
structure of random effects. The pair plot for the estimated 
random effects in the tree level is shown in Fig. 2. There 
was a weak positive correlation between the AGE and 
AGE random effects, but no substantial correlation 
between either of these random effects and the intercept 
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Fig. 2 Scatter plot of the estimated random effects in tree level from 
model 1.1 


random effects. A blocked diagonal matrix can be used to 
represent such covariance structure [12], with the intercept 
random effect corresponding to one block and the AGE and 
AGE random effects corresponding to another block. 
There was no remarkable correlation found among the 
random effects in the vertical stem position level. 

Several models with different structures for the vari¬ 
ance-covariance matrices of the estimated random effects 
were fitted and compared using the log-likelihood ratio test 
(LRT), the Akaike’s information criterion (AIC), and the 
Schwarz’s Bayesian information criterion (BIC). Accord¬ 
ing to the fit statistics presented on Table 2, model 1.2 with 
the blocked diagonal matrix in the tree level was the best of 
the variance-covariance structure of random effects. 
Hence, the variance-covariance structure could be repre¬ 
sented as 


Var (bi) = = 

000 

= 0 
0 


000 

0 

0 

0 

0n 

0 


0 

0n 

021 

0 

0 

022 


0 

012 , 

022 


Var (bij) = 'V 2 



Determining the structure of the within-group error 


The within-group error, e. hijk , were assumed to be inde¬ 
pendent for different /, j , k and independent of the random 
effects in model 1.2. The plots of residuals against the fitted 
values and other candidate variance covariates are useful 
for investigating within-group heteroscedasticity [12]. In 
this case, the cambial age is a natural candidate for the 
variance covariate. Figure 3 shows the plots of residuals 
versus AGE , indicating that the residuals decrease with 
AGE. Thus, we proceed by specifying the variance struc¬ 
ture of the within-group error to account for heterosced¬ 
asticity. We use a conditional error variance [12, 23], 
where we assume 

Var(fiy! | b u , ) = a 2 G 2 (n ijk , v ijk , d), (3) 


where is n ijk = E [y ijk I b Q , b t j], v ijk is a vector of variance 
covariates, 3 is a vector of variance parameters and G(.) is 
the variance function. A number of variance function can 
be used in the nlme package, the following two variance 
structures were tested in this study. The first is the power 
model which is given as 


Var(ep) = (j 2 \v ijk \ 25 , G( v ijk , <5) = 


Vijk 



The second is the exponential model which can be repre¬ 
sented as 


Var(ep) = a 2 exp(25vi jk ), G(v ijk , 5) = e\p(dv ijk ) (5) 

The parameter 3 is unrestricted, thus, the both variance 
structures can model a case where the variance increases or 
decreases with the variance covariate. There were signifi¬ 
cant increases in the log-restricted-likelihood, as evidenced 


Table 2 Comparisons of the model performance with different variance-covariance structures for the random effects 


Model 

Var-Cov structure 

Tree 

Position 

No. of parameters 

AIC 

BIC 

logLik 

LRT 

p value 

1.1 

Diag 

Diag 

13 

4808 

4850 

-2393.9 



1.2 

Block-diag 

Diag 

11 

4794 

4840 

-2386.1 

15.7 

0.00001 

1.3 

Diag 

Block-diag 

11 

4810 

4856 

-2393.9 

<0.001 

0.9996 

1.4 

Block-diag 

Block-diag 

12 

4796 

4864 

-2386.1 

15.7 

0.00004 


Diag diagonal matrix structure, Block-Diag blocked diagonal matrix structure, logLik log-restricted-likelihood estimated by restricted maximum 
likelihood method, AIC Akaike’s information criterion, BIC Schwarz’s Bayesian information criterion LRT likelihood ratio test calculated with 
respect to model 1.1 
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Fig. 3 Plot of residuals versus cambial age for the model 1.2 having 
homoscedastic within-group errors 

by the large value of the LRT, indicating that addition of 
the variance function to the model significantly improves 
on model 1.2 (Table 3). Based on AIC and BIC, power 
function of age (model 1.2.1) will be used to represent 
variance structure of the within-group error. 

Since the age trend of basic density could be considered 
as time series data, we need to pay attention to the temporal 
autocorrelation. Serial correlation structures are used to 
model dependence in time series data [24]. The empirical 
autocorrelation function provides a useful tool for investi¬ 
gating serial correlation in time series data. A plot of the 
estimated autocorrelation coefficients against lags for 
model 1.2.1 indicates that autocorrelations were significant 
even at cambial age lag 3 and 4 (data not shown). The 
following second-order moving average MA(2) was the 


best of the candidate correlation structures based on the 
statistics presented in Table 4. 

2 

o, = 52 Qj a t-j + ( 6 ) 

j =i 

where s t is the current within-subject error term, 0j are the 
moving-average parameters (j = 1, 2) and a t is a homo¬ 
scedastic noise term centered at 0 (E [a t ] = 0). The esti¬ 
mated normalized autocorrelation structure for model 
1.2.1.4 residuals was 

P = \f>{ 1 ), P{ 2 ), P{ 3 ), P( 4 ), P( 5 ), p(6),p(7),p(8)] T 

= [-0.198, -0.087, -0.112, -0.194, -0.098, 

-0.042, —0.097] T , 

( 7 ) 

where p(Z) is the empirical autocorrelation calculated at lag 
/ and the estimated parameters for the MA(2) model were 
= 0.348, and 9 2 = 0.256. 

We selected the MA(2) model that is more adequate to 
represent the autocorrelation at small lags even though the 
second-order autoregressive model AR(2) model and the 
MA(2) model indicated almost the same AIC, BIC and log- 
restricted-likelihood. This is because the empirical auto¬ 
correlation at larger lags tends to be less reliable due to the 
estimation by fewer residuals pairs [12]. In our data, there 
were differences of the ring number from pith among the 
disks. For example, a total ring number of a disk was 42 
while that of another was only 37, such disproportion of the 


Table 3 Comparisons of the model performance with different within-group error variance structures 


Model 

Variance function 

No. of parameters 

AIC 

BIC 

LogLik 

LRT 

p value 

1.2 

No structure 

11 

4794 

4840 

-2386.1 



1.2.1 

Power 

12 

4646 

4697 

-2311.2 

150.0 

<0.0001 

1.2.2 

Exponential 

12 

4682 

4732 

-2329.2 

113.8 

<0.0001 


Likelihood ratio test calculated with respect to model 1.2. Abbreviations are same in Table 2 


Table 4 Comparisons of the model performance with different within-groups correlation structures 


Model 

Correlation structure 

No. of parameters 

AIC 

BIC 

LogLik 

Test 

LRT 

p value 

1.2.1 

Independent 

12 

4646 

4697 

-2311.2 




1.2.1.1 

AR (1) 

13 

4623 

4678 

-2298.7 

1.2.1 versus 1.2.1.1 

25.0 

<0.0001 

1.2.1.2 

AR (2) 

14 

4621 

4680 

-2296.6 

1.2.1.1 versus 1.2.1.2 

4.1 

0.0419 

1.2.1.3 

MA (1) 

13 

4632 

4686 

-2303.0 

1.2.1.2 versus 1.2.1.3 

12.7 

0.0004 

1.2.1.4 

MA (2) 

14 

4621 

4680 

-2296.7 

1.2.1.3 versus 1.2.1.4 

12.6 

0.0004 

1.2.1.5 

ARMA (1,1) 

14 

4622 

4681 

-2297.2 




1.2.1.6 

ARMA (1,2) 

15 

4622 

4685 

-2296.0 

1.2.1.5 versus 1.2.1.6 

1.7 

0.1872 

1.2.1.7 

ARMA (2,1) 

15 

4623 

4685 

-2296.4 





AR autoregressive model, MA moving average correlation model, ARMA mixed autoregressive-moving average model, Other abbreviations are 
same in Table 2 
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Table 5 The estimated fixed effects parameters of basic density with 
the quadratic function of cambial age 


Parameters 

Estimate 

Std. error 

p value 

Po 

598.8215 

7.9846 

<0.0001 

Pi 

-13.1779 

0.7004 

<0.0001 

Pi 

0.1978 

0.0132 

<0.0001 


data also might be possibly the cause of the autocorrelation 
at larger lags. 

Final model 

We have tested the models having various mean and covari¬ 
ance structures, and have finally obtained the following model 
to describe the intra-tree variation in basic density. 


BDyk — (p 0 + hot + bo /j) + (/fi + bu + buj'j AGE\ 
+ (/?2 + ^2 i + ^2 ij) AGE\ -T Sijk, 




o 

o 

b 

0 

0 \ 

\ 



o 

o 

b 

0 

0 \ 

\ 

bj ~ N 

0 , 

0 


012 


, bij ~ N 

0 , 

0 

011 

0 




0 

021 

022 / 

/ 



0 

0 

022 / 

/ 


0 , 5) Hjjk((p, 8) G^(v m , 5)), 

G iJk (v ijk , 5) = \AGE\ S , H ijk (q>, 8) = ARMA(0, 2), 


( 8 ) 


where H ijk ((j), 6) is the serial correlation function and 
ARM A is mixed autoregressive-moving average model. 
The remaining elements of the model have been described 
previously. Parameter estimates, corresponding standard 
errors and p values for fixed effects of model 1.2.1.4 [Eq. 
(8)] are given in Table 5. 

A final assessment of the adequacy of model 1.2.1.4 [eq. 
(8)] is given by the plot of the augmented prediction [12, 
19] for a chosen tree (Fig. 4). Observed values of basic 
density, prediction by estimated fixed-effects parameters of 
model 1.2.1.4, which random effects were excluded, and 
that of containing random effects. The predictions con¬ 
taining the random effects follow the observed values 
closely, indicating that the final model explains the density 
variation data well. 

Discussion 

This study applied the linear mixed model to evaluate the 
effects of growth rate on the intra-tree variation in the basic 
density. The basic density decreased from pith to outward 
and the unique pattern was found in any height position 
(Fig. 1). The results were consistent with the previous 
reports in hinoki cypress [25, 26]. From Fig. 4, the final 
model with the quadratic function of cambial age [Eq. (8)] 













Cambial age 


Cambial age 


Cambial age 


Cambial age 


Fig. 4 Population prediction (setting all random effects to 0), 
subjects specified prediction (containing random effects) from the 
final model 1.2.1.4 and observed values of basic density versus 


cambial age. The solid lines, the dots line and filled circles indicate 
the population prediction, subjects specified prediction and observed 
values, respectively 
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explained the variation of basic density successfully. The 
density variation pattern found in hinoki is different from 
those of general coniferous tree, such as pine and larch, 
where the wood density increases with age [4, 7]. Decrease 
of wood density with age can be also found in sugi 
(Cryptomeria japonica) [26, 27]. 

There were no significant effects in the terms of plot, 
plot -AGE interaction, and plot -AGE interaction. The 
results indicate that growth rate does not affect the varia¬ 
tion pattern of basic density within the stem in hinoki 
cypress. There were few studies about the effects of growth 
rate on the variation pattern of basic density in hinoki 
cypress in spite of its importance for forest management 
and utilization. Fujiwara et al. [6] reported the similar 
results, but Hirai’s report [5] is inconsistent with our result. 
Growth rate affects wood density at different tree ages, and 
the inconsistency of these results would be due to the 
ignorance of age effects [4, 8]. The modeling approach 
employed in this study can assess the effects of growth rate 
on the variation of wood density quantitatively considering 
the age effects. The effect of growth rate on wood density 
varies greatly among species [4], and it can be confirmed 
that hinoki shows the similar tendency as the hard pines, 
Douglas-fir, and larch species. 

The defining characteristics of mixed-effects models 
are that they are applied to data where the observations 
are grouped according to one or more levels of experi¬ 
mental units and that they incorporate both fixed-effects 
terms and random-effects term [12]. Moreover, they 
present an inherent flexibility that allows for develop¬ 
ment of a unique variance-covariance structure alleviat¬ 
ing the problems of nonconstant variance and 
autocorrelation among the repeated measurements. The 
final model has two levels of mixed effects with random 
effects at the tree and vertical position levels. The ran¬ 
dom-effects estimates were found to be larger at the tree 
level than the vertical position level (data not shown). 
This means that the basic density values at each height 
position are relatively consistent, but the variation among 
trees is more noticeable, due to the unique patterns of 
basic density corresponding to height level. 

In this study, the variation of wood density within stem 
was modeled using the linear mixed-effects model. The 
methodology presented in this paper can easily extend for 
other wood properties [16, 18-20, 28]. In this case, it 
should be noted that the age trend of each traits is quite 
different. For instance, tensile strength in hinoki increases 
from pith to outward [29], and thus, the fitting model 
function should be different from that of wood density. We 
used the polynomial model that is linear in the parameters. 
By increasing the order of a polynomial model, one can get 
increasingly accurate approximations to the true regression 
function within the observed range of the data [12]. These 


empirical models are based only on the observed rela¬ 
tionship between the response and the covariates and do 
not include any theoretical considerations about the 
underlying mechanism producing the data. As the next 
step, it would be useful to apply the nonlinear model that is 
based on a model for the mechanism producing the 
response, and that also provides more reliable predictions 
for the response variable outside the observed range of the 
data [16, 19, 20]. 
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